A New Approach to Large-Scale Nuclear Structure Calculations 
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A new approach to large-scale nuclear structure calculations, based on the Density Matrix Renor- 
malization Group (DMRG) , is described. The method is tested in the context of a problem involving 
many identical nucleons constrained to move in a single large-j shell and interacting via a pairing 
plus quadrupole interaction. In cases in which exact diagonalization of the hamiltonian is possible, 
the method is able to reproduce the exact results for the ground state energy and the energies of 
low-lying excited states with extreme precision. Results are also presented for a model problem in 
which exact solution is not feasible. 
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PACS numbers: 21.60.Cs, 05.10.Cc 



The nuclear shell model Q is arguably the most pow- 
erful approach for a microscopic description of nuclear 
properties. In this approach, the low-energy structure of 
a given nucleus is described by assuming an inert doubly- 
magic core and then seeing how the effective interaction 
scatters the remaining nucleons over the valence orbits of 
the next major shell(s). Despite the enormous simplifica- 
tion provided by this shell-model approach, it is still only 
possible to describe nuclei in this way within limited re- 
gions of the periodic table, namely where the number of 
active nucleons or the degeneracy of the valence shells is 
sufficiently small. The most ambitious implementation 
of this method to date has been in a treatment of the 
binding energies of nuclei in the fp shell through ^"^Zn 

§. 

If we wish to use the shell-model approach in the de- 
scription of heavier nuclei or nuclei further from closed 
shells, we must come up with a reliable truncation pro- 
cedure, one that can reduce the number of shell-model 
configurations while maintaining the key dynamics of 
the interacting nucleons. Historically, many approaches 
have been used. Some truncate on the basis of weak- 
coupling considerations Q , others on the basis of symme- 
try considerations [Q , and others on the basis of Monte 
Carlo sampling Within the latter approach, it has 
recently proven possible to go beyond the fp shell to de- 
scribe the transition from spherical to deformed nuclei 
in the Barium isotopes In this work, we describe 
an alternative approach to large-scale nuclear structure 
calculations, called the Density Matrix Renormalization 
Group (DMRG). This method, developed originally in 
the framework of low-dimensional quantum lattice sys- 
tems 1^], was recently extended to finite Fermi systems 
The new methodology was first used in the treat- 
ment of a pairing problem Q of relevance to the physics 
of ultrasmall superconducting grains. Here we present 



its first application to a problem of relevance to nuclear 
structure. 

The basic idea of the DMRG method, as appropriate 
to finite Fermi systems, is to systematically take into ac- 
count the physics of all single-particle levels. This is done 
by first taking into account the most important levels, 
namely those that are nearest to the Fermi surface, and 
then gradually including the others in subsequent itera- 
tions. At each step of the procedure, truncation is imple- 
mented so as to optimally take into account the effects 
of the levels that are added while keeping the problem 
tractable. 

We will assume from the outset that each single- 
particle level in the problem admits the same number of 
possible states. This is the case, for example, when work- 
ing in a single-particle basis of axially-symmetric Nilsson- 
like levels, or if we consider pairs of time reversal states 
in a spherical basis. Each such level can accomodate four 
possible states, one with no particles, two with one par- 
ticle, and one with two particles. We denote the number 
of states that a given single-particle level admits as s; in 
this case, s = 4. 

Next, we assume that we have already treated some of 
the levels for particles and the same number for holes, 
namely those closest to the Fermi energy, and that the 
number of states in the two spaces is the same. We call 
that common number of states m. When we add the 
next particle level and the next hole level, the number of 
particle states increases from m to s x to and the number 
of hole states likewise increases from to to s x m. The 
DMRG method truncates from the s x to states to the 
optimum to states, both for particles and for holes. 

Following this optimal truncation, we then add the 
next levels for particles and holes and truncate again to 
the optimum to states for each. This procedure is contin- 
ued until all particle and hole levels have been sampled. 
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[Note: If one type of level is exhausted before the other, 
we subsequently add those levels that remain one at a 
time.] 

Finally, we carry out the calculation as a function of m, 
the number of particle and hole states that are kept. All 
previous applications of the DMRG method |^ have ex- 
hibited exponential convergence of the results (e.g., for 
the ground state energy) as a function of to, suggest- 
ing that this should likewise be the case in applications 
to nuclear structure. If so, we stop the procedure once 
the changes that arise with increasing to are acceptably 
small. 

The key question not yet addressed is "What do we 
mean by optimum and how do we implement a trunca- 
tion to those optimum states?" 

To answer these questions, we now consider the ground 
state of the full system, expressed as a sum of terms in- 
volving states in the particle space, \i >p, coupled to 
states in the hole space, \j >h , 

i=l,Np j=l,NH 

where Np is dimension of the particle subspace and Nh 
is the dimension of the hole subspace. [In the previous 
discussion, Np — Nh — s x to.] 

What we would like to do is to construct the optimal 
approximation to the ground state wave function l^* > 
that is achieved when we only retain to states in the par- 
ticle space and to states in the hole space. By optimal, 
we will mean that the projected wave function, the one 
that arises following the truncation, has the largest pos- 
sible overlap with the exact ground state wave function 
I* >. 

We will implement the truncation in two steps, first 
asking what is the best approximation when we truncate 
the particle states and then what is the best approxima- 
tion when we truncate the hole space. 

To arrive at the optimum truncation for particles, we 
first introduce the corresponding ground state reduced 
density matrix, 

3 = 1, Nh 

obtained by contracting over all the states of the hole 
space. We then diagonalize this Np x Np matrix, 

P^l^i" >p = u;P \u" >p . (2) 

A given eigenvalue oj^ represents the probability of find- 
ing the particle state |u" >p in the full ground state 
wave function of the system. The optimal truncation cor- 
responds to retaining the to eigenvectors that have the 
largest probability of being present in the ground state 
wave function, or equivalently those that correspond to 
the largest eigenvalues to^ ^ . 



Analogously, we construct the ground state reduced 
density matrix for holes, 

pfr- J2 (3) 

by contracting over particle states. If we diagonalize the 
Nh X Nh density matrix for holes and retain only the 
TO states with the largest eigenvalues, we are guaranteed 
to be choosing the best hole truncation in the sense of 
maximal overlap with the exact ground state. 

Summarizing, in each DMRG iteration we add to the 
system a new particle level and a new hole level. We 
then construct the haniiltonian matrix for the enlarged 
system and diagonalize it for the ground state and some 
low-lying excited states . From the ground-state wavefuc- 
tion, we calculate the reduced density matrices for parti- 
cles (|^) and holes (^) and diagonalize them. For each, we 
retain the to eigenvectors corresponding to the to largest 
eigenvalues. We then transform all operator matrices in 
the enlarged particle and hole spaces to the new trun- 
cated basis, completing the iteration. 

As noted earlier, the first application of this method- 
ology was reported in Ref. Q in the context of a pairing 
hamiltonian acting over a very large number of equally 
separated doubly-degenerate single-particle levels. The 
hamiltonian for this so-called picket fence model is 

^ = E '^J'^S^'^ ^3<^ ~ ^'^ S^- '^3'+ ' 

(4) 

with ej^ = j d and j — 1, 51 . Here, fl is the total 
number of levels and d is the spacing between adjacent 
levels The calculations of Ref. [^ assumed half filling, so 
that the total number of particles distributed over the 
levels is also il. 

As an example of the quality of the results that can be 
obtained with this method, consider the case of = 400. 
For this value of il, the dimension of the full hamilto- 
nian matrix that would have to be diagonalized is 10^^^ ! 
While this problem is obviously much too large to be 
solved by standard diagonalization techniques, it can 
be solved to arbitrary accuracy using a method devel- 
oped by Richardson in the 1960s [|l^. For a problem in 
which A — 0.224, the Richardson solution for the ground 
state has a correlation energy of -22.5183141 in miits of 
d. When the same problem is solved using the DMRG 
method with to = 60, a ground state correlation energy 
of -22.5168 in the same units is achieved. The agreement 
is to better than 1 part in 10**, despite the fact that the 
maximum dimension hamiltonian matrix that had to be 
treated for this value of to was only 3066. 

Clearly, for a pairing hamiltonian acting in a uniform 
doubly-degenerate single-particle space, the Density Ma- 
trix Renormalization Group method works remarkably 
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well. On the other hand, the hamiltonian (Q) of this 
problem is extremely simple, being equivalent to a one- 
body hamiltonian for a hard-core boson. In that respect, 
even though the previous results are suggestive that the 
method might also work well for problems in nuclear 
structure physics, we still need to demonstrate this more 
convincingly by applying it to a fermion problem with a 
true two-body interaction. 

With that in mind, we have now completed the first 
test of the new DMRG methodology in nuclear structure. 
We considered a schematic model in which a large num- 
ber of identical particles are restricted to a single large-j 
shell and interact via a sum of a pairing plus quadrupole 
force. Since a single j shell does not give rise naturally 
to a Fermi surface, we also included a single-particle en- 
ergy term in the hamiltonian, one that favors an oblate 
solution. The hamiltonian of the model is 

H = -xQ-Q- gP^ P - e ^ |m| c]„c,„ . (5) 

m 

Because of the last term, the hamiltonian is not rotation- 
ally invariant, so that its eigenstates do not have good 
angular momentum. 

The first results we present are for 10 particles in a 
j ~ 25/2 orbit. The dimension of the hamiltonian ma- 
trix that would have to be diagonalized (in the m scheme) 
for this problem is 109,583, which can be readily handled 
using the Lanczos algorithm. 

In Table I, we present the results for x = Ij <7 = 0, 
and e = 0.1. The first row gives the exact ground state 
energy; subsequent rows give the ground state energy 
obtained using the DMRG approach as a function of m. 
At the end of each row, we show the maximum dimen- 
sion hamiltonian matrix that must be diagonalized in the 
DMRG procedure. 

As is evident from the table, the results converge very 
rapidly to the exact ground state energy. By m—QQ, we 
obtained a result that is off by only 1 part in 10^. For 
this choice of m, the largest matrix we had to diagonalize 
was 398 X 398. 

To see what happens in the presence of both 
quadrupole and pairing correlations, we next present the 
results for X = 1, 9 — 0.05, and e = 0.1. Such a value for 
the pairing strength leads to a strong depletion of prob- 
ability from the hole levels, but still leaves a well defined 
Fermi surface. 

Results for the ground state energy of the system are 
summarized in Table II. The exact result is —16.00367, 
obtained using the Lanczos algorithm. In the DMRG ap- 
proach, we get an energy of —16.00151 for m=80, with 
a maximum hamiltonian dimension of 685. The results 
continue to improve slightly with increasing m, as the 
maximum dimension of the hamiltonian matrix increases. 

The DMRG results for the ground state energy are 
plotted in Fig. 1 as a function of m. Included in the 



figure is an exponential fit to these results. When ex- 
trapolated, this exponential predicts an asymptotic end 
result of —16. 00358 ± .00027, in excellent agreement with 
the exact ground state energy. 

Table III shows results for the excitation energies of 
the lowest three excited states. The agreement for the 
excited states is almost as good as for the ground state, 
even though the method, as described earlier, only tar- 
geted the ground state in the optimization procedure. 

The bottom line is that for these calculations, in which 
we could compare with the results of Lanczos diagonal- 
ization, we obtain excellent agreement with the exact re- 
sults, not just for the ground state but for higher ex- 
cited states as well. Furthermore, the excellent results 
are achieved while diagonalizing matrices of moderate di- 
mensions. 

The fact that we could achieve a high level of accuracy 
not just for the ground state but for low-lying excited 
states as well may be a reflection that all these states 
have the same intrinsic structure. When different intrin- 
sic structures enter, it may be necessary to modify our 
optimization criterion to include mixed density matrices 
that contain information on more than one state of the 
system [g. 

The excellent quality of the results we obtained for a 
j = 25/2 orbit encouraged us to treat a more complex 
system, one for which exact diagonalization is not possi- 
ble. We considered the case of a j = 55/2 orbit with 20 
particles. The other parameters of the calculation were 
X = 1, g = 0.1, and e = 0.1. In this case, the exact cal- 
culation would involve a matrix of dimension 5.31064 x 
10^^. The results for the ground state energy are shown 
in Table IV. By m = 60, the calculations have clearly 
converged and we should have a reliable ground state 
energy to about six significant figures. 

In our view, these calculations demonstrate very 
clearly the great promise of the DMRG method in nu- 
clear structure. As long as there is a well defined Fermi 
surface in the problem, the method leads to extremely 
accurate results not only for the ground state but for 
low-lying excited states as well. Furthermore, the fairly 
rapid convergence that typically arises as a function of m 
suggests that the method can be used quite reliably for 
very large-scale calculations. 

It is worthwhile here to expand briefly on the rationale 
for the assumption that all levels of the single-particle 
basis admit the same number of states. This assumption 
greatly facilitates implementation of an iterative scheme 
in which the addition of new levels can be readily acco- 
modated with no change of formalism. Relaxation of this 
assumption may be feasible, but no doubt at a significant 
cost to computational simplicity. 

Now that we have completed this first test of the 
methodology with such impressive success, we are plan- 
ning to gradually expand the complexity of the problems 
we consider. Our ultimate goal of course is to use this 
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method to treat very large-scale nuclear structure prob- 
lems involving both neutrons and protons populating a 
set of non-degenerate single-particle levels and interact- 
ing via a general interaction. 
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Table 1: Ground state energy for 10 particles in a j=25/2 level. 
The Hamiltonian parameters are: x = 1, g = 0,and e = 0.1. 



m 




Egs 




Max{Dim) 


Exact 




-15.58837 




109,583 


20 




-15.58798 




106 


40 




-15.58830 




217 


60 




-15.58836 




398 


80 




-15.58836 




656 


Table 2: Ground state energy for 10 particles in a j= 


=25/2 level. 




The Hamiltonian parameters are: x 


= 1,9 = 0.05, 


e= .1. 




m 




Egs 




Max{Dim) 


Exact 




-16.00367 




109,583 


20 




-15.99380 




118 


40 




-15.99776 




238 


60 




-15.99946 




548 


80 




-16.00151 




685 


100 




-16.00246 




964 


120 




-16.00271 




1267 


140 




-16.00306 




1648 


Table 3: Excitation energies for 10 particles 


in a J = 25/2 level. 


The hamiltonian parameters are: 


X = l,g = 0.05 




and e = 0.1. 








m 


El 




E2 




Exact 


0.51643 




0.87141 


1.10294 


40 


0.52464 




0.96894 


1.18556 


60 


0.51774 




0.90385 


1.11528 


80 


0.51854 




0.88581 


1.11353 


100 


0.51817 




0.88274 


1.11432 


120 


0.51772 




0.88068 


1.11115 


140 


0.51743 




0.87842 


1.10942 


Table 4: Ground state energy for 20 particles in a j= 


=55/2 level. 




The Hamiltonian parameters are: x 


= 1,5 = 0.1, 


e= .1. 




m 




Egs 




Max{Dim) 


Exact 




? 




5.31X 10''' 


20 




-103.98844 




100 


30 




-103.99420 




180 


40 




-103.99574 




240 


50 




-103.99827 




361 


60 




-103.99894 




430 
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